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ABSTRACT 

We use three-dimensional smoothed particle hydrodynamics simulations to study the 
non-linear thin shell instability (NTSI) in supersonic colliding flows. We show that for 
flows with monochromatic perturbations and for flows with white-noise perturbations, 
growth speeds approximate quite well to the analytic predictions of Vishniacl ( [T99 4). 



For flows with subsonic turbulence, growth speeds match Vishniac s predictions 
only at short wavelengths where the turbulence is weaker. We find that supersonic 
turbulence, of a lower Mach number than the colliding flows, completely suppresses 
the NTSI. Our results provide a diagnostic for identifying the presence of the NTSI in 
colliding flows with turbulence. 
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1 INTRODUCTION 

Colliding flows of gas are common in the interstellar medium 
and lead to the formation of dense, shock-compressed layers 
(e.g. [Draine & McKee]|1993 i. Since stars form from dense 
gas, colliding flows may be important triggers of star formation. 
Colliding flows can occur at the confluence of expanding H II re- 



gions (^Elmegreen & Lada|1977j , stellar-wind bubbles l,Stevens, 
Blondin & Pollock||1992) ~and supernova remn ants ([Williams 



et al.||1997|, or from large-scale galactic flows ffionnell et al 



2om 



Vazquez-Semadeni et al. 



Burkert||2008| [Hennebelle et al. 



2006[ |Heitsch, Hartmann & 
2008 1. Turbulence, which is 



believed to exist at a range of scales in the Galaxy (e.g. |Mac| 
|Low & Klessen|2004[ l, also creates colliding flows. 

In this paper we use three-dimensional smoothed particle 
hydrodynamics (SPH) simulations to explore the effect of 
the non-linear thin shell instability (NTSI; |Vishniac]|1994| l on 
the development of a shock-compressed layer between two 
supersonically colliding flows. We seed the inflowing gas with 
a variety of initial perturbations, and compare the resulting 
NTSI growth speeds with the analytic predictions of Vishniac 
|T994|. Previous studies have been limited to two-dimensional 
simulations ffilondin & Marks|| 1996| [Klein & W6odsl[T998l 



ondin & Marks 



attempted to 



[Hueckst aedt 20031, and only 
measure the dependence of growth speeds on wavenumber. 



1.1 NTSI 

The NTSI occurs in a dense layer that is already corrugated, 
as illustrated in Fig. [T] The thermal pressure exerted by the 
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Figure 1. The non-linear thin shell instability. Gas flows in from the 
left and right, confining the layer with ram pressure (large arrows). The 
dense gas in the layer resists the compression with thermal pressure 
(small arrows), but this is always normal to the surface of the layer. 
The misalignment of ram and thermal pressure causes shear in the layer 
(shown with arrows within the layer), leading to density enhancements 
(shown in dark grey). 



material in the layer always acts perpendicular to the surface 
of the layer, whereas the ram pressure of the inflowing gas 
acts along the collision axis. The resulting shear drives material 
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towards the extremities of the corrugation, and this causes the 
perturbation to grow. 

The NTSI was first described by [Vishniac] ( |1994^ , who 
considered colliding flows in an isothermal gas. He derived an 
expression for the growth time T as a function of the wavelength 
A, the amplitude rj and the sound speed Cj, 



-r(A,r7) 



(1) 



It follows that short-wavelength modes grow fastest, and the 
growth speed is given by 



T(A,ri) 



J7N3/2 



(2) 



The NTSI is a non-linear instability and can only grow if 
it is seeded by a finite perturbation with an initial amplitude T]o 
comparable to or larger than the thickness of the layer, 

^(0^ — , (3) 
where vq is the speed of the inflowing gas. 



1.1.1 Range of unstable wavelengths 

[Vishniac] jl994^ predicts that the NTSI can grow for only a 
limited range of wavelengths. For a stationary slab in a highly 
supersonic collision, the large wavelength limit is given by 



(4) 



which states that wavelengths are not excited on larger scales 
than those coordinated by sound waves. The small wavelength 
limit is given by the thickness of the layer; a thick layer cannot 
be effectively corrugated by a small wavelength perturbation, 
and so 



1.1.2 Growth to saturation 



(5) 



Once the NTSI has been excited at a wavelength X, it grows to 
saturation at a growth speed jVishniac|I994^ 



3/2 



<n<ic. 



(6) 



The lower limit is the growth speed of a perturbation with an 
amplitude equal to the layer thickness, which is the smallest 
amplitude perturbation unstable to the NTSI. The higher limit 
can be rearranged to give a growth time limit, 

X 



(7) 



which states that the NTSI does not grow in less than the sound 
crossing time and is another expression of the large wavelength 
limit of equation |4] Equation |6] also implies that the growth 
speed is not strongly supersonic; perturbations where r] ^ A 
will saturate due to the large bending angle. 

Initially, smaller wavelength perturbations grow fastest. As 
the layer thickness grows, smaller wavelength perturbation stop 
growing and are overtaken by larger perturbations. The fastest- 
growing mode has a wavelength of approximately the layer 
thickness, and will subsequently saturate due to the increasing 
thickness of the layer. 



Table 1. Simulation parameters 



Cross-section 








^end 


'«SPH 


A'SPH 




(pc) 


(pc) 


(pc) 


(Myr) 


(Mq) 




Narrow 


1 


0.5 


0.01 


0.1 


9.23 X lO-' 


8 X 10'' 


Square 


2 


0.5 


0.5 


0.15 


1.48 X lO-* 


5 X 10« 



1.2 Simulation parameters 

We collide flows of isothermal molecular gas of temperature 
lOK and sound speed 0.19kms^'. The gas has an initially 
uniform density of lO^^^gcm^^. The flows collide along the 
X-axis with a collision velocity vq of ±3.8kms^', and so Mach 
number of 20. All simulations are three-dimensional, and 
are periodic in the y and z directions. 

The collision box has dimensions L^, Ly and along 
the X-, y- and z-axes, respectively. We conduct two types of 
simulation: simulations with a narrow rectangular cross-section 
and one-dimensional perturbations, and simulations with a 
square cross-section and either two-dimensional perturbations 
or turbulence. The narrow cross-section simulations are of a 
higher resolution than the square cross-section simulations. 
Table [T] lists the box size for these simulations, together with 
the simulation end time /endi mass per SPH particle nispH ™d 
total number of SPH particles A'sph- 

We impose one-dimensional perturbations which vary 
across the v-axis, and two-dimensional perturbations which 
vary across the y- and z-axes. In this work we use dimensionless 
wavenumbers measured in terms of the relevant box size, which 
is always 0.5 pc. Wavenumbers are therefore given by 



0.5 pc 



(8) 



These dimensionless wavenumbers can be converted to physical 
wavenumbers by multiplying by 2pc^' . 

We impose three types of initial velocity perturbation on 
these flows: monochromatic sinusoidal perturbations, white- 
noise perturbations, and an initial turbulent velocity field. The 
white-noise perturbations are formed from a superposition of 
sinusoidal perturbations of equal amplitude and random phase. 
Monochromatic and white-noise perturbations are only imposed 
on a small region close to the contact discontinuity between the 
colliding flows, reducing to zero at a distance ^peit away from 
the contact discontinuity by a smoothing function 





+ 1 




V "pert / 





if \x\ < R, 



pert 



(9) 



,0 



if Ul > R 



pert 



The region containing the perturbation, and its corresponding 
momentum, is rapidly accreted onto the layer, after which only 
unperturbed gas is accreted onto the layer. The instability is 
then free to grow and decay without interference from other 
perturbations. 

For turbulent perturbations the initial turbulent velocity 
field covers the entire simulation. The layer experiences a rela- 
tively constant input of perturbations as the turbulent medium is 
accreted onto the layer. Any growth of the NTSI must compete 
with these incoming perturbations. 
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1.2.1 Monochromatic perturbations 

For a single monochromatic sinusoidal perturbation of integer 
wavenumber k and amplitude A, the additional velocity due to 
the perturbation is given by 



(10) 



0.2 



8vx{x,y,z) = vqA sin ( ^rikj- ) f{x) , 



for a one-dimensional perturbation. In two dimensions we use a 
sinusoidal square-grid perturbation given by 



Svx{x,y,z) = vqA sin (^nkyj- 



sin Ink, 



(11) 



where we have integer wavenumbers ky and k^; for monochro- 
matic perturbations we use ky — k^. This scjuare-grid pertur- 
bation corresponds to a superposition of the four plane-wave 
modes k = {±ky,±k,) with equal amplitude (but varying sign); 
we consider only the k = (ky,k^) mode in our analysis. As 
described in Section [2.3.3| we have 



B + k^ . 



(12) 



and so when ky = the wavenumber ky is related to the 
wavenumber khy k = \/2ky. 

This means that it is not possible to excite the same 
wavenumbers of one-dimensional and two-dimensional pertur- 
bations, as the periodic boundary conditions constrain us to use 
integer wavenumbers of k for one-dimensional perturbations, 
and integer wavenumbers of ky and k~ for two-dimensional 
perturbations. For the same set of integers, the wavenumbers for 
the simulations with two-dimensional perturbations are larger, 
and the wavelengths correspondingly smaller 

We use an amplitude A of 0.2 (except for simulations with 
two-dimensional monochromatic perturbations where we use an 
amplitude of 0.5), and a perturbation extent i^pert (see equa- 
tion |9j of 0.05 pc. We perform simulations using monochro- 
matic wavenumbers k or ky of 2, 4, 6, 8, 10, 12, 14 and 
16 for one-dimensional and two-dimensional monochromatic 
perturbations. Figure |2] shows an example of a simulation with 
a one-dimensional monochromatic perturbation. 



1.2.2 White-noise perturbations 

White-noise perturbations are generated by the superposition of 
sinusoidal perturbations between k = I and k = 20 for simula- 
tions with one-dimensional perturbations, and all combinations 
of ky and k, between 1 and 20 for simulations with two- 
dimensional perturbations. Each wavenumber is given the same 
amplitude A, and a random phase of between and 27t chosen 
from a uniform distribution. We repeat the process to generate 
ten random sets of perturbations. Figure[3]shows an example of 
a simulation with one-dimensional white-noise perturbations. 



1.2.3 Turbulent perturbations 

In our work we use only divergence-free turbulence generated 
from a Fourier spectrum with a power-law form A(k) oc k^^ 
(or equivalently a power spectrum P{k) o= k^'^). We generate 
the cubic and periodic velocity field in a similar way to that 
described by |Mac Low et al.| (|1998). The turbulent velocity 
field is imposed at the start of the simulation and subsequently 
decays. Turbulent fields are only applied to the lower-resolution 



4x10'^° 
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Figure 2. Den,sity cross-section plot of simulation with one-dimensional 
monochromatic k=12 perturbation after 0.09 Myr. Gas is inflowing from 
the left and right, driving the NTSI in the dense shock-compressed layer. 



simulations with square cross-sections, where they are scaled to 
fit across the y and z dimensions, and are then repeated across 
the X dimension as required. This preserves the periodic nature 
of the y- and z-axes. 

Ten velocity-field realizations are generated. We scale each 
velocity field to two different average velocities to give a 
subsonic velocity field with Mach number 0.5 and a supersonic 
velocity field of Mach number 4. This gives ten subsonic and ten 
supersonic velocity fields. In both cases the turbulence is only 
a small perturbation on the colliding flows, which have Mach 
number 20. 

Figure|4]shows an example of a simulation with supersonic 
turbulence. The layer appears uniformly displaced to x > 
in this figure; this is due to large-scale perturbations in the 
turbulence. Slices at different z coordinates show displacement 
or bending of a similar magnitude but both left and right of 
x = 0. 



2 METHOD 



We use the SEREN ( [Hubber et al.|2011[ l code, which solves the 
equations of fluid dynamics using the SPH method l |Gingold"&l 
Mona ghan|1977[[Lucy|197'7jl. This cod e includes a conservative 
'grad-h' SPH scheme ('Springel & Hernquist 2002j [Price &| 
|Monaghan||2004| l. We do not include self-gravity or magnetic 
fields. 



2.1 Creation of initial conditions 

For both the narrow and square cross-section simulations, we 
generate five realizations of particle positions. We randomly 
place particles in a periodic box which covers the desired sizes 
Ly and L, in the y- and z-axes, respectively, to ensure there 
are no regularly repeating patterns of SPH particles along these 
axes. However, to save computational time, the periodic box has 
only a fraction of the desired size in the x-axis. The particles 
are allowed to settle under the influence of hydrodynamics and 
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Figure 3. Density cross-section plot of simulation with one-dimensional 
white-noise perturbations after 0.05 Myr. Gas is inflowing from the left 
and right. 

artificial viscosity. Following the settling procedure, the particle 
positions in the box are repeated along the x-axis to generate the 
desired box length L^. 

For the simulations with one-dimensional monochromatic 
perturbations, each of the monochromatic wavenumbers 2, 4, 6, 
8, 10, 12, 14 and 16 is combined with each of the five settled 
distributions to create a total of forty sets of initial conditions. 
This is repeated for the simulations with two-dimensional 
monochromatic perturbations. For the white-noise simulations, 
there are ten random sets of velocity perturbations; each of the 
five settled particle distributions are therefore used twice. The 
same applies for the ten subsonic and ten supersonic simulations 
with turbulence. 



2.2 Gridding of particle data 

The SPH data is interpolated on to a three-dimensional cubic 
grid of density using the SPLASH package ( |Price|[2007] l. This 
grid is of size 256"^, is centred at the centre of the collision 
box, and has sides of length 0.5 pc. Some gas is therefore not 
included in the grid, and in the simulations with narrow cross- 
sections the grid covers some empty space across the z-axis. 

For narrow cross-section simulations, which contain only 
one-dimensional perturbations, we calculate the average gas 
density in the simulation box across the z-axis. This reduces the 
three-dimensional grid to a two-dimensional grid of size 256^. 

Each simulation produces a number of snapshot outputs at 
regular intervals of 0.005 Myr. The gridding process is repeated 
for each output snapshot to build a time-sequence of density 
grids for each simulation. 

2.3 Identification of layer 

Once the data has been gridded, cells containing the dense 
shock-compressed layer are identified. A density threshold of 
5 X 10^^' gcm^"* is used to identify potential 'layer' cells. In 
simulations containing supersonic turbulence, dense regions can 



form which are distant from the layer; to remove these from our 
list of layer cells, we first have to estimate some layer properties. 

We calculate the median x position of layer cells along 
each column of grid cells running along the collision axis. This 
gives an estimate for the current position of the layer. For each 
column, we calculate the median absolute deviation of layer 
cells X position from this estimated layer position; this gives 
an estimated layer half-thickness for the column. Finally, we 
take the median of the estimated layer half-thickness for each 
column to produce an estimated half-thickness for the entire 
layer. 

We then identify discrete groups of connected poten- 
tial layer cells. We find the average position of each group 
{Gx,Gy,Gz). We then find the estimated x position of the layer 
aty = Gy and for square cross-section simulations z = G~. The 
difference between this position and Gj gives an estimate of the 
perpendicular distance of the group from the layer G±, which 
we compare with the size of the group along the x-axis, Gl. 

If Gx is larger than Gl by more than the estimated layer 
half-thickness, the group is rejected. This is approximately the 
same as saying that the group is rejected if it does not overlap 
the layer. All cells within the rejected group are removed from 
our list of potential layer cells. We find this simple filtering 
routine, which accounts for variations in the layer position 
but assumes a constant layer thickness, works well to remove 
unrelated dense regions formed by turbulence without affecting 
the identification of the layer. 

2. 3. 1 Layer position 

Having filtered the layer, we now calculate more accurately the 
layer position, including only 'layer' cells in our analysis. The 
layer position for each grid column along the x-axis is given 
by the centre-of-mass along that column. This produces a one- 
dimensional grid of layer positions for narrow cross-section 
simulations and a two-dimensional grid of layer positions for 
square cross-section simulations. We use these layer positions 
as our tracer of bending-mode instabilities present in the layer. 

2.3.2 Layer thickness 

In order to calculate the maximum resolvable wavenumber as a 
function of time, we need to calculate or estimate the thickness 
of the layer. We define the layer thickness for each column in 
X as the total extent of layer cells along x for that column. This 
includes any empty voids in the layer, as we find this is a better 
representation of the true thickness of the layer This is later 
adjusted to account for the bending angle of the layer. 

2.3.3 Fourier transform 

To obtain the Fourier spectrum of perturbations in the layer, 
we take the Fourier transform of the grid of layer positions, 
yielding a one-dimensional Fourier spectrum for the narrow 
cross-section simulations and a two-dimensional grid of Fourier 
amplitudes for the square cross-section simulations. This is 
normalized such that values represent the amplitude of pertur- 
bations in parsecs. 

For the square cross-section simulations, we perform 
a simple circular averaging technique to convert the two- 
dimensional grid of Fourier amplitudes, defined at integer 
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Figure 4. Density cross-section plot of simulation with supersonic 
turbulence after 0.07 Myr. Gas is inflowing from the left and right. 



wavenumbers fc,, and k~, into a one-dimensional Fourier 
spectrum. This is defined over a similar range of wavenumbers 
to the Fourier spectrum for the narrow cross-section 
simulations, but is defined at a large number of fractional 
wavenumbers rather than only at integer wavenumbers. A 
perturbation with wavenumbers and will create a peak at 

the one-dimensional wavenumber k= \/^? + k}- 



2.3.4 Calculation of growth speeds 

For each simulation, we have calculated a time-sequence of 
one-dimensional Fourier spectra. We use the central differential 
method to calculate the growth speed of perturbations in layer 
position in units of pcMyr^'. This yields the growth speed as 
a function of time for each wavenumber. As we have multiple 
realizations of each simulation, we average growth speeds 
across these realizations, and calculate the standard deviation. 

For monochromatic simulations, we are only interested in 
the growth speed for the wavenumber which was originally ex- 
cited, averaged across the five realizations per wavenumber. We 
therefore extract the growth speed as a function of time for this 
single wavenumber. For square cross-section monochromatic 
simulations, we can directly extract the relevant wavenumber 
from the two-dimensional Fourier transform without using 
our circular averaging. We combine the results from all the 
monochromatic simulations at different wavenumbers to pro- 
duce the growth speed as a function of time and as a function of 
all the excited wavenumbers. 



3 RESULTS 

The collision of the flows creates a shock-compressed layer. As 
expected, the thickness of this layer increases linearly with time, 
and the rate of increase does not depend strongly on the type 
of initial velocity perturbation. However, the rate does differ 
between narrow and square cross-section simulations. We fit 
the thickness as a function of time for each simulation using 




0.04 0.06 
Time (Myr) 



0.10 



Figure 5. The growth speeds for one-dimensional monochromatic 
perturbations as a function of time. Each line represents a different 
wavenumber; error bars show the standard deviation across realizations. 



the linear function i^fit(^) = ™d average over the different 
simulations to find K = 0.47pcMyr^' for narrow cross-section 
simulations and K = 0.22pcMyr^' for square cross-section 
simulations. This is faster than the 0.019pcMyr^' predicted for 
an idealized isothermal collision. 

This increase is not due to insufficient resolution, as 
the layer grows faster for the higher-resolution narrow cross- 
section simulations than the lower-resolution square cross- 
section simulations. It is due to the bending of the layer, which 
reduces the effective ram pressure. This causes the layer to reach 
a lower density than for the idealized isothermal case. 

For a layer at an angle 6, the component of the inflowing 
velocity vo that contributes to the ram pressure is vq cos 0. Since 
ram pressure is given by pv^, the effective ram pressure and 
therefore the layer density is reduced by a factor of cos^ 9. This 
increases the thickness of the layer by a factor of sec^ 6. Both 
our initial perturbation and the NTSI create significant bending 
angles, which greatly increase the thickness of the layer. 

This effect has been identified previously in simulations 
of the NTSI. |Blondin & Marks] ( |1996| l found their layers grew 
faster than expected even in the absence of the NTSI due to 
the growth of the Kelvin-Helmholtz instability. In simulations 
of cold dense layers confined by supersonic shocks, |Folini &] 
|Walder| ^2006^ found that all their layers reached a maximum 
density contrast relative to the inflowing medium of '-^ 30 
regardless of the Mach number. This is similar to our layers, 
which have average density contrasts of approximately 16 
and 35 for the narrow and square cross-section simulations, 
respectively. 

Figure |5] shows the growth speed of different modes 
as a function of time for simulations with one-dimensional 
monochromatic perturbations. The results for each mode are 
taken as the average from a set of simulations where only 
that mode is perturbed. As expected, the growth speeds of all 
wavenumbers begin at a value corresponding to the amplitude 
of the initial perturbation, and then quickly decay. 

Following the decay of the initial perturbation, the growth 
speeds of perturbations at small wavenumbers do not increase 
appreciably during the simulation time. However, higher 
wavenumber perturbations start to grow, and reach larger 
growth speeds than the initial perturbation, before decreasing 



© 2013 RAS, MNRAS QOQ.pTfn] 



6 A. D. McLeod and A. P. Whitworth 



again at later times. This increase is due to tlie NTSI; the 
later decline is due to the layer becoming too thick to support 
perturbations of higher wavenumbers. 

We can compare the growth of the NTSI to the maxi- 
mum and minimum wavelength limits defined by Equations [4] 
and|5] respectively, and the corresponding minimum wavenum- 
ber limit 



1 



and the maximum resolvable wavenumber 



(13) 



(14) 



We find the minimum wavenumber limit is not reliable. 
At the end of the simulation time the minimum dimensionless 
wavenumber should be much greater than 26, and higher at 
earlier times. However, Figure|5]suggests a minimum wavenum- 
ber closer to 10. The minimum wavenumber limit assumes that 
wavenumbers do not grow in less than the sound crossing time 
across a wavelength. As described in Section [l.l.2[ combining 
this with the limit of small bending angles also implies that 
the growth speeds of perturbations should remain at most 
weakly supersonic. We note that several of these modes achieve 
growth speeds which imply supersonic velocities at the peaks of 
the perturbation, but none reach more than approximately Mach 
2.5. 

However, we do confirm the maximum wavenumber limit, 
equation[T4] We use the averaged rate of layer thickness growth 
described above, with a different value for narrow and square 
cross-section simulations, to estimate the layer thickness and 
the corresponding maximum resolvable wavenumber at each 
time. In our simulations, growth at each wavenumber stops and 
decays close to the predicted value. 



3.1 Monochromatic perturbations 

For monochromatic perturbations, we first consider simula- 
tions with a one-dimensional perturbation where we excite 
a single wavenumber using the initial velocity perturbation 
given by equation [TO] We then consider simulations with a 
two-dimensional perturbation, where we excite a single set of 
wavenumbers using the velocity perturbation given by equa- 
tion[TT| As noted in Section [T.2.1| we cannot excite exactly the 
same wavenumbers in these simulations. 



3.1.1 One-dimensional perturbations 

Figure |6] shows the growth speeds for simulations with one- 
dimensional monochromatic perturbations. As expected, at 
early times the slope of the line is roughly flat, as all modes are 
given the same initial perturbation amplitude. At later times, 
higher wavenumbers grow faster than smaller wavenumbers, 
appearing to ap proximately match the slope predicted by 



Vishniac 



(1994} of A:' -". The highest wavenumbers begin to 
decay at the latest times, but at somewhat higher wavenumbers 
than predicted by equation[T4| 

Figure [7] shows the index of power-law fits to the growth 
speeds as a function of wavenumber. We fit from the smallest 
wavenumbers to the maximum resolvable wavenumber at each 
time. As expected, the index begins at approximately zero, but 
rises rapidly. This rise brings it somewhat above the predicted 
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Figure 6. The growth speeds for one-dimensional monochromatic 
perturbations as a function of wavenumber. Each line represents results 
from a different timestep; the upper frame displays the growth speed 
at early times and the lower frame displays the growth speeds at 
later times. Error bars show the standard deviation across realizations. 
The gradient of the solid black line indicates the relation predicted 
by |Vishniac| (1994) . The an'ows indicate the maximum resolvable 
wavenumber for the correspondingly coloured timestep; arrows run 
from latest to earliest timestep going left to right, and some arrows may 
be at higher wavenumbers than can be plotted. 



index of 1.5, instead lying in the range of approximately 2.0 to 
2.5. 



3.1.2 Two-dimensional perturbations 

Figure [8] shows the growth speeds for simulations with two- 
dimensional monochromatic perturbations. Again the slopes 
begin fairly flat before rising to approximately match the predic- 
tion of|Vishniac| At later times the decay of higher wavenum- 
bers is evident and occurs at approximately the wavenumber 
predicted by equation [T?) 

Figure |9] shows the index of power-law fits to the growth 
speeds as a function of wavenumber. This initially rises from 
zero, as expected, and then remains close to the predicted value 
of 1.5. 
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Figure 7. Indices of power-law fits to growth speeds of one-dimensional 
monochromatic perturbations as a function of time. The number of 
points available to construct each fit is also shown; fits constructed 
with fewer points are less reliable. The black dashed line indicates the 
prediction of |Vishniacljl994| . 
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Figure 9. Indices of power-law fits to growth speeds from simulations 
with two-dimensional monochromatic perturbations, as for Figure|7] 

3.2 White-noise perturbations 

White-noise perturbations are created by a superposition of 
velocity perturbations of equal amplitude and random phase, as 
described in Section [T.2.2| Not all modes can grow in any given 
simulation; some modes will grow while others will decay, 
creating the large standard deviations shown in Figures [T0| 
and[T2] 

Unlike for monochromatic perturbations, where we con- 
sider only a single wavenumber from each simulation, we 
consider the growth speed of a range of wavenumbers and 
wavevectors from each simulation. For simulations with two- 
dimensional perturbations, this means we have to perform 
circular averaging, as described in Section [2.3.3[ to produce the 
growth speed as a function of the one-dimensional wavenumber. 
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Figure 8. The growth speeds from simulations with two-dimensional 
monochromatic perturbations as a function of wavenumber, as for 
Figure|6] 



3.2.1 One-dimensional perturbations 

Figure [To] shows the growth speeds for simulations with one- 
dimensional white-noise perturbations. Again the slopes begin 
fairly flat before rising to approximately match the prediction 
of |Vishniac| At later times, the decay of higher wavenumbers 
is evident and occurs at or slightly above the wavenumber 
predicted by equation[T4] 

Figure [TT] shows the index of power-law fits to the growth 
speeds as a function of wavenumber. As expected this initially 
rises from zero and then remains close to the predicted value of 
1.5, although with considerable scatter The results at later times 
should be ignored due to the low number of wavenumbers that 
can still be resolved as the layer grows in thickness. 

3.2.2 Two-dimensional perturbations 

Figure [T2] shows the growth speeds for simulations with two- 
dimensional white-noise perturbations. Due to the circular av- 
eraging, results are presented as a continuous line with a shaded 
error zone instead of individual points. These growth speeds 
are more complex than for the simulations with monochromatic 
perturbations. It is clear that higher wavenumbers generally 
grow faster up to some limiting wavenumber However, there 
is a large decrease in growth speeds at smaller wavenumbers 
at early times, followed by a growth at small wavenumbers at 
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Figure 10. The growth speeds from simulations with one-dimensional 
white-noise perturbations as a function of wavenumber, as for Figure|6] 
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Figure 12. The growth speeds from simulations with two-dimensional 
white-noise perturbations as a function of wavenumber, as for Figure|6] 
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Figure 11. Indices of power-law fits to growth speeds from simulations 
with one-dimensional white-noise perturbations, as for Figure|7] 




0.02 



0.04 0.06 0.08 

Time (Myr) 



Figure 13. Indices of power-law fits to growth speeds from simulations 
with two-dimensional white-noise perturbations, as for Figure|7] 



later times. Growth speeds appear to decay at wavenumbers 
somewhat smaller than predicted by equation[T4| 

Figure [T3] shows the index of power-law fits to the growth 
speeds as a function of wavenumber. This initially rises from 
zero as expected, before reaching approximately 2. The index 
then decreases with time, falling to approximately 0.7. Our 



power-law fits may not accurately reflect the growth of the NTSI 
for two reasons: the unexpected growth at small wavenumbers 
and the decay at larger wavenumbers below the maximum 
wavenumber predicted by equation[T4] 
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Figure 14. The growth speeds for simulations with supersonic turbu- 
lence as a function of wavenumber, as for Figure |6] 



3.3 Turbulence 

We now consider simulations that do not have an imposed 
sinusoidal velocity perturbation, but instead have a turbulent 
velocity field imposed on all the gas in the simulation. We use 
only square cross-section simulations, with either subsonic or 
supersonic turbulence as described in Section [1.2.3| Although 
the turbulence itself is three-dimensional, it provides a time- 
varying two-dimensional set of perturbations as it is accreted 
onto the layer. 

When fitting power-law slopes to the growth speeds, we 
fit from a minimum wavenumber of = 8 to help separate the 
effect of the turbulence and the NTSI. We fit to the maximum 
resolvable wavenumber predicted by equation [14] or = 64 if 
this is smaller. 



3.3.1 Supersonic turbulence 

Figure [14] shows the growth speeds for simulations with su- 
personic turbulence. The influence of the turbulence, which is 
strongest at small wavenumbers, dominates the growth speeds. 
No evidence for the NTSI can be seen, with growth speeds 
decreasing with increasing wavenumber. 

Figure [15] shows the index of power-law fits to the growth 
speeds as a function of wavenumber. The index begins strongly 
negative, in contrast to earlier simulations where the index 
begins at approximately zero. This is because there is not equal 
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Figure 15. Indices of power-law fits to growth speeds from simulations 
with supersonic turbulence, as for Figure|7] 
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Figure 16. The growth speeds for simulations with subsonic turbulence 
as a function of wavenumber, as for Figure|6] 



amplitude at all modes; the turbulence is stronger at smaller 
wavenumbers. The index remains at approximately —2. There 
is thus no evidence for the NTSI. 
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wavenumbers. They found the dependence on wavenumber 




shghtly weaker than predicted, finding T 



, while we 
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Figure 17. Indices of power-law fits to growth speeds from simulations 
with subsonic turbulence, as for Figure|7] 



3.3.2 Subsonic turbulence 

Figure [T6] shows the growth speeds for simulations with sub- 
sonic turbulence. The turbulence remains dominant at smaller 
wavenumbers; however, the growth speeds are smaller by an 
order of magnitude. At higher wavenumbers a rise can be seen 
not too dissimilar in slope to the predictions of |Vishniac|jl994^ . 
Qualitatively there is a range of wavenumbers over which 
growth speeds increase with increasing wavenumber until the 
limit imposed by the thickness of the layer. 

Figure [TT] shows the index of power-law fits to the growth 
speeds as a function of wavenumber. Due to the complex 
nature of these growth speeds, with turbulence dominating at 
small wavenumbers and the NTSI dominating at intermediate 
wavenumbers, such fits are difficult. As previously described, 
for our turbulent simulations we ignore wavenumbers below 
/: = 8 in our fitting. This excludes the range of wavenum- 
bers still dominated by the turbulence, where growth speeds 
decrease with increasing wavenumber, and includes only the 
region where growth speeds are increasing with increasing 
wavenumber up to the maximum resolvable wavenumber. 

As for the supersonic turbulence, the initial index is ap- 
proximately —2. However, this rapidly rises, eventually peaking 
around 0.7. Although this does not quantitatively match the 
predicted index of |Vishniac|jl994| >, it is a qualitative match in 
that higher wavenumbers grow more quickly. 



find the dependence slightly stronger than predicted. Qualita- 
tively our results are in good agreement; the NTSI grows expo- 
nentially until saturation, which occurs when the layer thickness 
exceeds half the perturbation wavelength. The layer remains 
turbulent and bloated, and is much thicker than expected for 
an idealized isothermal shock, similar to the results of [Folini &] 
|Walder|p006| >. 

[Blondin & Marks] began their simulations with a sinu- 
soidally perturbed slab already in place, and did not use an 
initial velocity perturbation. However, their simulations, using 
a similar Mach number of 25, are otherwise similar to our 
simulations with a one-dimensional monochromatic perturba- 
tion. Our results therefore extend their work to two-dimensional 
perturbations, white noise and turbulent initial conditions, but 
we find their basic results hold for all cases we examine except 
those containing turbulence. 

I Klein & Woods] ([T998 l studied the growth of the NTSI in 



colliding clouds, while Hueckstaedt 



(|2003|l conducted a study 



of colliding flows similar to that of Blondin & Marks] jT996| 
Both used two-dimensional axisymmetric codes. As in the 
simulations of [Hunter et al.| ( |1986 ), both found that the NTSI 
requires strong cooling. When using an adiabatic equation of 
state, ,Klein & Woods| found the NTSI was not excited; they 
were also able to show that the NTSI did not grow in the ab- 
sence of initial perturbations, as predicted by | Vishniac| ^1994^ 
Both sets of simulations are in qualitative agreement with our 
results; however, neither made a quantitative measurement of 
the dependence of growth rates on wavenumber. 



4.1 NTSI in the Galaxy 

It is worth considering what astrophysical systems could be 
susceptible to the NTSI. The NTSI requires a layer with a sig- 
nificant density contrast, bounded by shocks from inflowing gas 
on both sides. The inflowing gas must therefore be supersonic 
and must cool rapidly in the shock. We also find that supersonic 
turbulence can prevent the NTSI from being observed if the 
growth speeds of turbulent perturbations accreted onto the 
layer exceed those of the NTSI. Since supersonic turbulence is 
ubiquitous in the interstellar medium, the NTSI may be replaced 
in real systems by these turbulent perturbations, which would 
limit the applicability of our results. 



4 DISCUSSION 



jHunter et al.| jI986^ simulated the collision of supersonic 
flows using a two-dimensional grid code. They found that the 
dense layer formed became unstable to a 'Rayleigh-Taylor- 
like' instability, and noted that it only occurred in the presence 
of cooling. This instability, the NTSI, was first analysed for 
supersonic isothermal collisions by | Vishnlac] { 1 994) . The NTSI 
requires a large density contrast between the inflowing and layer 
gas, and so requires strong cooling to produce an approximately 

isothermal shock, as found by |Hunter et al.|(T986). 

Using a two-dimensional Eulerian code, Blondin & Marks| 
( |1996| l studied the growth of the NTSI in a layer formed 
between planar flows. These simulations confirmed the analysis 
of [Vishniac] ( |1994[ ) that growth rates were larger at larger 



4.1.1 Colliding clouds 

Colliding molecular clouds are likely to provide the conditions 
for the NTSI, provided the collisions are not excessively off- 
centre, and are likely to be the systems to which our simulations 
are most applicable. The numerical study of colliding clouds 
has a long history, beginning with the one-dimensional and 
two-dimensional sim ulations of ] Stone|(|1970a|b^. Many fiirther 
studies followed (e.g. |Smith|I980||Hausman|1981[|Gilden|1984[ 



Latta nzio et al.|I985) , however many of these would be unable 
to identify the NTSI due to low resolution, the use of a mirror 
boundary instead of a second cloud, or the lack of initial 
perturbations. The exception is the work by jKlein & Woods| 
( |1998| l described earlier. 

Simulations including cooling showed that such colli- 
sions are approximately isothermal (e.g. |Stone|I970a[|Hausmaril 
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|1981[ |Sabano & Tosa|[l985| l. Together with the supersonic 
motions of clouds, this provides the high density contrast 
required for growth of the NTSI. More recent simulations such 
as those by |Kitsionas & Whitworth| ( |2007 |l have tended to focus 
on the resulting star formation rather than instabilities that occur 
during the collision. 

Both subsonic and supersonic turbulence are common in 
molecular clouds and cores. Previous studies have not consid- 
ered the effect of turbulence on the NTSI; our results for sub- 
sonic and supersonic turbulence can be applied to the collisions 
of such turbulent clouds. Further numerical work studying the 
quantitative growth of the NTSI in colliding clouds is needed. 



4.1.2 Large-scale flows 

Converging flows on galactic scales have been suggested as the 
precursors to molecular clouds, since such flows can collect 
material in a relatively dense layer jHeitsch et al.|2005[ |2006[ 
[Vazquez-Semadeni et al . 2006 2007|. Such flows may trigger 
the NTSI and the Kelvin-Helmholtz instability in the layer, 
which drives turbulence in the layer. 

However, unlike in our simulations the inflowing gas is hot 
atomic gas. When this gas is compressed it undergoes a thermal 
instability to form a cold, dense layer. This means that the 
dense layer and any fragments formed are confined by thermal 
pressure as well as ram pressure, and so our results may not 
be directly applicable. Nonetheless the NTSI is an important 
process if the collision is sufficiently supersonic and the cooling 
is su fficiently rapid jHeitsch et al.|2005) . 

[Audit & Hennebelle (2010i simulated converging flows 
using either a two-phase medium with a thermal instability or an 
isothermal equation of state. The velocity of the inflowing gas 
is the same in both cases, but corresponds to a Mach number 
of 1.5 for the two-phase medium and of the order of 10 for 
the isothermal case. The two-phase medium leads to a much 
more fragmentary structure, with many small objects which are 
compressed by the surrounding hot medium. The isothermal 
case is very similar to our simulations of the NTSI, with a more 
clearly defined layer and many bending modes. 



4.1.3 Colliding winds 

Converging flows occur in the colliding wind binaries examined 
by [Stevens et al.| jl992( l. They find that the shocked layer 
produced in such collisions can be isothermal for binaries with 
a period of less than a few days, particularly for the denser 
winds of Wolf-Rayet stars. For adiabatic collisions the shocked 
layer is unstable to the Kelvin-Helmholtz instability, but when 
the collision is isothermal then a thin shell instability disrupts 
the layer; |Vishniac| ( |1994| l identified this as the NTSI. The 
physical conditions of these collisions is quite different to those 
of our simulations, but provided the collision is sufficiently 
radiative the dependence of growth rates on wavenumber should 
be similar in the region where the collision is close to planar. 



4.1.4 Other converging flows 

The NTSI may also occur in the bow shocks of runaway stars, 
but our results are less likely to be applicable here; even if the 
bow shocks are isothermal the dominant instability is likely 
to be the transverse acceleration instability and not the NTSI 



JDgani, van Buren & Noriega-Crespo 1996a bl. Similarly, the 
confluence of expanding H II regions ( Elmegreen & Lada| I977] l 
or supernovae remnants ( [Williams et al.| 19*97 1 also provide the 
converging flows required for the NTSI, but may not provide the 
strong radiative cooling needed for the NTSI. 



4.2 Observations 

We can speculate as to possible observational evidence of the 
NTSI in real systems. Our results are most closely applicable 
to molecular cloud collisions and so we focus on this case. 
We study the growth speeds of perturbations by identifying the 
position of the dense layer between the clouds. While this may 
be possible for clouds that collide exactly parallel to the plane 
of the sky, for most systems it may be easier to observe the 
line-of-sight velocity of layer perturbations using molecular line 
observations. 

For collisions along our line of sight, the dense layer can 
be observed using a suitable optically-thin molecular tracer. If 
the clouds are sufficiently supersonic, their emission should 
not overlap in velocity space with the dense, mostly subsonic 
layer. If some of cloud emission is at similar velocities to the 
layer, then a tracer with a critical density higher than the cloud 
density but lower than the layer density could be used. This 
is complicated by the relatively low density contrast between 
layer and inflowing gas we observe; another possible alternative 
would be a tracer of shocked gas such as SiO. 

Provided suitably high-resolution observations are avail- 
able then the velocity of the layer provides equivalent informa- 
tion to the layer displacement for our purposes. The velocity 
across the collision plane can be processed in the same way as 
we process the layer displacement across the collision plane, 
except no time derivative is required. 

We analyse the growth speeds of perturbations by com- 
puting the Fourier spectrum of layer displacements and taking 
the time derivative of the Fourier spectra; this is equivalent to 
taking the time derivative of the layer displacements, i.e. the 
layer velocity, and computing the Fourier spectrum. This will 
require observations of the velocity across the collision plane; 
as an example, we estimate the required resolution for a 10 pc 
cloud at varying distances. 

For a minimal Fourier spectrum, we can require a linear 
resolution of 50 points across the collision plane, which will 
provide information to a wavenumber of 25. There is likely to be 
a density gradient across the dense layer from centre to edge due 
to the spherical nature of the clouds, so we restrict observations 
to the central lOpc width and so require a linear resolution of 
0.2 pc. 

The Atacama Large Millimeter/submillimeter Arrray 
(ALMA), in its widest configuration, has a maximum spatial 
resolution of between 6 and 37 mas. Assuming a resolution of 
20 mas, this yields our desired linear resolution out to a distance 
of several megaparsec, and so we easily achieve the required 
spatial resolution for any target in the Milky Way. More 
difficult is the required spectral resolution. Since our maximum 
velocity is of the order of the sound speed, approximately 
0.2 km s^', we require a velocity resolution much smaller than 
this. The highest spectral resolution available from ALMA is 
approximately O.Olkms^', which only gives a total velocity 
resolution of about 20 points over the subsonic regime. 

We can speculate as to the time that the NTSI could be 
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observed in such a collision. Once the NTSI is saturated, its 
effects can be seen as a bloated, turbulent layer, but our results 
will not be applicable to this phase of the layer evolution. 
In order to observe the slope of the Fourier spectrum, we 
require smaller wavelength modes to still be growing, which 
requires that the layer is thinner than half the wavelength. For 
a wavenumber of 20 across the central 10 pc width, the layer 
must be thinner than 0.25 pc. Assuming a Mach number of 20, 
and a density contrast of 30, this will last for approximately 
1 Myr. This is similar to the collision time for the clouds, so for 
this case the NTSI may be growing during the majority of the 
collision. 

4.3 Magnetic fields 

In this work we do not explore the effect of magnetic fields. 
Magnetic pressure can support a layer and reduce the density 
contrast between a layer and the inflowing medium (|Stone| 
|1970a[ ISabano & Tosa|[l985| l, which should decrease the ef- 
fectiveness of the NTSI. The magnetized NTSI was studied in 
detail by [Heitsch et al.| ( |2007^ for Mach 4 flows, who found 
magnetic fields normally weaken the NTSI either by resisting 
shear, if the field is aligned with the inflow, or by stiffening 
the equation of state, for fields perpendicular to the inflow. 
The effect of magnetic fields will depend on the balance of 
ram pressure and magnetic field strength, and thus may be less 
important for our large Mach number flows. 



5 CONCLUSIONS 

Our simulations have explored the effect of initial perturbations 
on the growth of the NTSI. We have defined a method to 
identify the position of a dense shock-compressed layer within a 
simulation, identify the perturbations within that layer, calculate 
the growth speeds and so help identify the presence of the NTSI. 
We have found that, in general, one-dimensional perturbations 
produce a similar result to two-dimensional perturbations, sug- 
gesting that the two-dimensional NTSI is similar to the one- 
dimensional NTSI analysed by | Vishnlac] ( | 1 994] ) . 

We have found that the thickness of the layer increases 
faster than for the idealized case, but does not depend much 
on the type of the initial perturbation. This is not a simple 
problem of resolution, since the layer thickness increases faster 
in the higher-resolution narrow cross-section simulations than 
in the much lower-resolution square cross-section simulations. 
The thickness of the layer does provide a good estimate of the 
maximum resolvable wavenumber as a function of time, usually 
to within a factor of two. 

We have only considered a single initial collision velocity 
and single initial density, and therefore have the same range 
of wavenumbers susceptible to the NTSI in each case. A more 
complete study might examine the effect of these parameters on 
the minimum wavenumber susceptible to the NTSI. We have 
found the original minimum wavenumber limit proposed by 
|Vislmiac] ( |1994^ to be unreliable, as the layer perturbations can 
achieve both supersonic speed and high bending angles. 

For simulations with one-dimensional monochromatic per- 
turbations, we have qualitatively matched the predictions of 
[Vishniac] (| 1994*1, although the slope of the growth speeds is 
slightly higher than predicted. For two-dimensional monochro- 
matic perturbations, however, the slope of the growth speeds 



is consistent with 'Vishniac's analytic prediction. It should 
be noted that this prediction was limited to one-dimensional 
monochromatic perturbations. 

Simulations with one-dimensional white-noise perturba- 
tions also match the analytic prediction. Simulations with two- 
dimensional white-noise perturbations are more complex, but 
match the general pattern that higher wavenumbers grow more 
rapidly. 

We also consider the effect of turbulence, which is likely to 
be found in realistic colliding flows. We find that for supersonic 
turbulence, the turbulence dominates over the NTSI, and the 
growth speeds reflect the Fourier amplitude spectrum of the 
turbulence. The growth speeds for the simulations with subsonic 
turbulence are similar at small wavenumbers, except an order of 
magnitude weaker. This allows the growth speed of the NTSI to 
exceed the growth speed of turbulent perturbations over a range 
of larger wavenumbers, as the turbulence is weaker at larger 
wavenumbers. 

Our results provide a diagnostic to identify the growth 
of the NTSI in the presence of turbulence. In the absence 
of the NTSI, growth speeds will decrease with increasing 
wavenumber if, as is normally assumed, the turbulence is 
stronger at smaller wavenumbers. If the NTSI is present, there 
will be a range of wavenumbers over which growth speeds will 
increase with increasing wavenumber. The limits of this range 
are the smallest wavenumber where the NTSI grows faster than 
turbulent perturbations and the largest wavenumber where the 
NTSI can grow (set by the thickness of the layer). 

This diagnostic may be used to identify the presence of 
the NTSI in colliding flows in more complex and more realis- 
tic simulations of astrophysical phenomena, such as colliding 
molecular clouds. It may also be of use to identify the NTSI in 
observations of molecular gas in the interstellar medium, where 
turbulence is believed to be ubiquitous. 



ONLINE MATERIALS 

Movies of a representative sample of our simulations can be 
found online accompanying this paper on the Monthly Notices 
website. For one-dimensional perturbations, there is a simula- 
tion with a = 8 monochromatic perturbation, and a simulation 
with white noise. For two-dimensional perturbations, there is 
a simulation with a two-dimensional ky = 8 monochromatic 
perturbation, and a simulation with two-dimensional white 
noise perturbations. There is also a simulation with subsonic 
turbulence and a simulation with supersonic turbulence. There 
are density cross-sections for all simulations; for simulations 
with two-dimensional perturbations or turbulence, column den- 
sity along the x-axis is included in separate movie files. 



ACKNOWLEDGMENTS 

Plots of SPH data in this work were produced with SPLASH 
( |Price|[2007^ . This research has made use of NASA's Astro- 
physics Data System. Simulations were performed using the 
Cardiff University ARCCA (Advanced Research Computing 
CArdiff) cluster Merlin. Andrew McLeod was funded during 
this work by an STFC Studentship at Cardiff University. The 
authors would like to thank the anonymous referee for their 
helpful and constructive comments. 



© 2013 RAS, MNRAS 000,[T|fT3l 



Simulations of the non-linear thin shell instability 13 



REFERENCES 

Audit E., Hennebelle P., 2010, A&A, 511, A76 

Blondin J. M., Marks B. S., 1996, New A, 1, 235 

Bonnell I. A., Dobbs C. L., Robitaille T. P., Pringle J. E., 2006, 

MNRAS, 365, 37 
Dgani R., van Buren D., Noriega-Crespo A., 1996a, ApJ, 461, 372 
Dgani R., van Buren D., Noriega-Crespo A., 1996b, ApJ, 461, 927 
Draine B. T., McKee C. R, 1993, ARA&A, 31, 373 
Elmegreen B. G., Lada C. J., 1977, ApJ, 214, 725 
Folini D., Walder R., 2006, A&A, 459, 1 
Gilden D. L., 1984, ApJ, 279, 335 

Gingold R. A., Monaghan J. J., 1977, MNRAS, 181, 375 
Hausman M. A., 1981, ApJ, 245, 72 

Heitsch R, Burkert A., Hartmann L. W., Slyz A. D., Devriendt J. E. G., 

2005, ApJ, 633, LI 13 

Heitsch R, Slyz A. D., Devriendt J. E. G., Hartmann L. W., Burkert 

A., 2006, ApJ, 648, 1052 
Heitsch R, Slyz A. D., Devriendt J. E. G., Hartmann L. W., Burkert 

A., 2007, ApJ, 665, 445 
Heitsch R, Hartmann L. W., Burkert A., 2008, ApJ, 683, 786 
Hennebelle P., Banerjee R., Vazquez-Semadeni E., Klessen R. S., 

Audit E., 2008, A&A, 486, L43 
Hubber D. A., Batty C. R, McLeod A., Whitworth A. P, 201 1, A&A, 

529, A27 

Hueckstaedt R. M., 2003, New A, 8, 295 

Hunter J. J. H., Sandford I. M. T., Whitaker R. W., Klein R. I., 1986, 
ApJ, 305, 309 

Kitsionas S., Whitworth A. P, 2007, MNRAS, 378, 507 
Klein R. I., Woods D. T., 1998, ApJ, 497, 777 

Lattanzio J. C. Monaghan J. J., Pongracic H., Schwarz M. P., 1985, 

MNRAS, 215, 125 
Lucy L. B., 1977, AJ, 82, 1013 

Mac Low M., Klessen R. S., Burkert A., Smith M. D., 1998, Physical 

Review Letters, 80, 2754 
Mac Low M.-M., Klessen R. S., 2004, Reviews of Modern Physics, 

76, 125 

Price D. J., 2007, PASA, 24, 159 

Price D. J., Monaghan J. J., 2004, MNRAS, 348, 139 

Sabano Y., Tosa M., 1985, Ap&SS, 1 15, 85 

Smith J., 1980, ApJ, 238, 842 

Springel V., Hernquist L., 2002, MNRAS, 333, 649 

Stevens I. R., Blondin J. M., Pollock A. M. T., 1992, ApJ, 386, 265 

Stone M. E., 1970a, ApJ, 159, 277 

Stone M. E., 1970b, ApJ, 159, 293 

Vazquez-Semadeni E., Ryu D., Passot T., Gonzalez R. F., Gazol A., 

2006, ApJ, 643, 245 

Vazquez-Semadeni E., Gomez G. C, Jappsen A. K., Ballesteros- 

Paredes J., Gonzalez R. R, Klessen R. S., 2007, ApJ, 657, 870 
Vishniac E. T., 1994, ApJ. 428, 186 

Williams R. M., Chu Y.-H., Dickel J. R., Beyer R., Petre R., Smith 
R. C, Milne D. K., 1997, ApJ, 480, 618 

This paper has been typeset from a T^X/ MgX file prepared by 
the author. 



© 2013 RAS, MNRAS Q0Q,pTfT3] 



